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Abstract 

On the occasion of its 25th anniversary, we revise Chen’s derivative rule for electron tunneling [C. 
J. Chen, Phys. Rev. B 42, 8841 (1990)] for the purpose of computationally efficient simulations of 
scanning tunneling microscopy (STM) based on first principles electronic structure data. The revised 
model allows the weighting of tunneling matrix elements of different tip orbital characters by an 
arbitrary energy independent choice or based on energy dependent weighting coefficients obtained 
by an expansion of the tip single electron wavefunctions/density of states projected onto the tip 
apex atom. Tip-orbital interference in the STM junction is included in the model by construction 
and can be analyzed quantitatively. As a further advantage, arbitrary tip geometrical orientations 
are included in the revised model by rotating the coordinate system of the tip apex using Euler 
angles and redefining the weighting coefficients of the tunneling matrix elements. We demonstrate 
the reliability of the model by applying it to two functionalized surfaces of recent interest where 
quantum interference effects play an important role in the STM imaging process: N-doped graphene 
and a magnetic M 112 H complex on the Ag(lll) surface. We find that the proposed tunneling model 
is 25 times faster than the Bardeen method concerning computational time, while maintaining good 
agreement. Our results show that the electronic structure of the tip has a considerable effect on STM 
images, and the Tersoff-Hamann model does not always provide sufficient results in view of quantum 
interference effects. For both studied surfaces we highlight the importance of interference between 
s and p z tip orbitals that can cause a significant contrast change in the STM images. Our method, 
thus, provides a fast and reliable tool for calculating STM images based on Chen’s derivative rule 
taking into account the electronic structure and local geometry of the tip apex. 

PACS numbers: 68.37.Ef, 71.15.-m, 73.63.-b 


1 


I. INTRODUCTION 


The role of the tip in scanning tunneling microscopy (STM) on the imaging contrast mech¬ 
anisms has been extensively studied in recent years. A reliable interpretation of experimental 
findings can be obtained by STM simulations that need electron tunneling models capable of 
dealing with a diversity of tip parameters in a consistent manner. Furthermore, it is required 
that STM simulations are computationally inexpensive and user-friendly in order to provide 
a useful tool not only for theoretician experts of electronic structure methods but also for 
experimental STM groups. 


The pioneering electron tunneling model was proposed by Bardeen who derived the tunnel 


a 


ing current formula based on first order perturbation theory |1|. While this method contains 
the effect of the electronic structure of the STM tip explicitly, it is only applicable in the 
pure tunneling regime and does not account for multiple scattering effects. An extension of 
the Bardeen method including multiple scattering has been proposed by Palotas and Hofer 


that has been implemented in the BSKAN code [3|. Using the Keldysh nonequilibrium 
Green’s function formalism, Mingo et al. developed an electron transport model that is valid 
in both tunneling and close-contact regimes [4|. These methods require the single electron 
wavefunctions or Green’s functions of the sample surface and the tip or of the coupled sys¬ 
tem that can be obtained from first principles electronic structure calculations. Although 
these single electron quantities can routinely be calculated, e.g., by using density functional 
theory (DFT), STM image simulations are still considered as computationally demanding if 
performed at high levels of electron transport theory. The reason is the numerous parameters 
affecting the electron tunneling that have to be considered if aiming at an accurate interpreta¬ 
tion of STM experiments, e.g., the bias voltage, the tunneling current or tip-sample distance 
and the tip characteristics: material, geometry, functionalization, termination, orientation, 
and the corresponding electronic structure. Note that the modeling of the tip needs special 
effort as the above listed tip parameters can be chosen in practically infinite combinations. 

To overcome the computational drawback of high level electron transport theories, several 


simplifications have been introduced. Tersoff and Hamann assumed an s-wave tip flfl. This 
way the Bardeen current formula has been recast to contain the electronic structure of the 
sample surface only. In the Tersoff-Hamann model the tunneling current is proportional to 
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the integrated local density of states of the sample at the position of the tip. The energy 
integral corresponds to the bias voltage. This approach has become the most widely used 
method for simulating STM and scanning tunneling spectroscopy (STS) due to its relative 
simplicity, although it completely neglects the effects of the STM tip. 


Allowing more generality for the tip orbitals, Chen derived the so-called derivative rule 
from the Bardeen current formula, where the resulting tunneling matrix element is propor¬ 
tional to the linear combination of spatial derivatives of single electron wavefunctions of the 


sample surface in the vacuum [7]. The spatial derivatives are determined by the orbital char¬ 
acters of the STM tip. However, practical calculations of the coefficients for the linear combi¬ 
nation were not reported. The determination of the energy dependent coefficients in Chen’s 
derivative rule from first principles calculations is one of the subjects of the present work. 
Palotas et al. developed a tunneling model using the three-dimensional Wentzel-Kramers- 
Brillouin (3D-WKB) theory |8j] that accounts for the energy dependent combination of tip 
orbital characters determined from the orbitally decomposed DOS projected to the tip apex 
atom 9). We will show that the orbital-decomposed projected tip DOS can be related to an 
approximation of the tunneling matrix weighting coefficients in the revised Chen’s derivative 
rule. 


The above mentioned tunneling models have been successfully used in theoretical inves¬ 
tigations of a multitude of STM junctions. To select a few examples, we focus on studies 
concerned with the role of the tip geometry, orbital character and functionalization on the 
STM imaging. Teobaldi et al. investigated the effect of bias voltage and tip electronic 
structure on the STM contrast formation of the highly oriented pyrolytic graphite (HOPG) 
surface using tungsten tips with different terminations and sharpnesses jlO]. Chaika et al. 
demonstrated that by using oriented single crystalline tungsten tips it is possible to select 
a particular tip electron orbital for high-resolution imaging of HOPG [111. Il2j . Channel 
selective tunneling was also examined by Wong et al. using tip functionalization with hexa- 
peri-hexabenzocoronene (HBC) molecules [13]. Employing Chen’s derivative rule, Gross et al. 
simulated STM images of pentacene and naphthalocyanine molecules using CO-functionalized 
tips m| . The increased lateral resolution achieved by these tips demonstrated the significant 
contribution of p-type tip states. Siegert et al. studied the influence of s- and p-wave tip 
symmetries on the STM maps of 7T-conjugated molecules using the reduced density matrix 
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formalism combined with Chen’s derivative rule 15]. The orbital-dependent 3D-WKB model 
was extended for spin-polarized STM in Ref. jl6|. The effect of tip orbital symmetries on the 
scanning tunneling spectra was also investigated by probing the cuprate high-temperature 
superconductor Bi2Sr2CaCu208+<5. Suominen et al. found that the symmetry of the tip can 
radically change the topographic image due to the overlap of sample and tip orbitals |l7| . 


while da Silva Neto et al. stated that the apparent nematic behavior of the lattice is likely 
related to a realistic STM tip probing the band structure of the material [18]. They also 


pointed out the importance of tunneling interference effects in the STM junction. 

The effect of inter- and intra-atomic interference of electron orbitals has also been in the 
focus of several studies. Telychko et al. investigated the nitrogen-doped graphene surface with 
tungsten and diamond tips and found significantly smaller current ('dip') above the nitrogen 
atom than above the neighboring carbon atoms at constant height [19]. This finding has been 
explained by a destructive quantum interference essentially resulting from the C-N 7T-bond. 
Using the Keldysh Green’s function formalism, Jurczyszyn and Stankiewicz and Mingo et al. 
extensively investigated inter-orbital interference effects in various tip-sample combinations 
and found that the interference has a considerable influence on STM images and STS spectra 

20, 211. Sachse et al. showed that an antiferromagnetic alignment of Mn spin moments in 


a M 112 H complex on the Ag(lll) surface explains the experimental STM observation of a dip 


above the middle of the Mn dimer 


22], In the present work we point out that a destructive 


quantum interference between s and p z tip orbitals contributes to the emergence of such a 
dip in the STM image. 

Very recently, the role of the spatial orientation of the tip in the tunneling process has been 
the subject of a few studies. Hagelaar et al. performed STM simulations of NO adsorbed on 
Rli(lll) surface over a wide range of tip apex terminations and orientations and compared 


them with experimental STM images 


23|. They found that asymmetric tip orientations 


provide good qualitative agreement with the experiments under certain tunneling conditions. 
Lakin et al. developed a technique to recover the relative orientation of a Cgo-functionalized 
tip and a Ceo molecule adsorbed on the Si(lll)-(7x7) surface based on Chen’s derivative 
rule I 24 ]]. Using the 3D-WKB method, Mandi et al. studied tip rotational effects on the 
imaging of HOPG(OOOl) 25] and W(110) surfaces 26] and found that the STM images can 
considerably be distorted due to different spatial orientations of the tip. Although the 3D- 
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WKB method lacks the inclusion of electron interference, it proved to be useful for statistically 
determining the geometry and likely orientations of the tip in bias voltage dependent STM 
experiments of HOPG based on previously inaccessible large scale simulations of tip models 
in a computationally efficient way [27|. The conclusion of all these listed studies is that 


the relative orientation of the sample surface and the local tip apex geometry is far from 
being highly symmetric that is usually considered in standard STM simulations nowadays. 
Therefore, it is clear that spatial orientations of the tip have to be integrated into STM 
simulation models, thus, the second subject of the present work is to include arbitrary tip 
geometrical orientations into Chen’s derivative rule. 

The paper is organized as follows. The revised Chen’s derivative rule is presented in section 
im Computational details of the electronic structure calculations of surface and tip models 
are reported in section IIIII Results of simulated STM images and comparisons with other 
electron tunneling methods are found in section HV Al and IIV Bl for the N-doped graphene and 
for the antiferromagnetic Mn 2 ff complex on Ag(lll) surface, respectively. We summarize 
our findings in section [V] 


II. REVISED CHEN’S DERIVATIVE RULE 

In the tunneling regime the current I depending on the bias voltage V between the sample 
surface and the tip can be calculated using Bardeen’s tunneling formula jl|, 

HV) = [1 - /(£„ + e V)} \M^\ 2 5(E„ -E„- eV), (1) 

fll/ 

where e is the elementary charge, h is the reduced Planck constant, / is the Fermi distribution 
function and M^ is the tunneling matrix element between two single electron states of the 
sample (/i) and the tip ( v ) involved in the tunneling. E M and E v denote corresponding Kohn- 
Sham eigenenergies that can be obtained from first principles calculations. Note that /i and 
v denote composite indices of the band (n), wave vector (ky) and spin (a) in the separate 
sample and tip subsystems, respectively. The tunneling matrix element can be calculated as 
an integral over the S separation surface in the vacuum between the sample and the tip, 

M,e = -%-( dS. (2) 

J S 
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In the model the tunneling is assumed to be elastic and energy conservation is ensured by the 
Dirac-delta in Eq. (JT]). At finite temperature the thermal broadening of the electron states 
has to be taken into account. This is usually done by approximating the Dirac-delta with a 
Gaussian function, 


S{E p - E v 


eV) ~ 


1 

= exp 


{E„ — E v — eVf 
2A 2 


( 3 ) 


In principle, all /i — u transitions have to be considered with the probability given by this 
Gaussian factor (and \M pu \ 2 ) when calculating the tunneling current but practically transi¬ 
tions with significantly low probability can be neglected, e.g., if E p — E v — eV\ > 3A, where 
A = ksT is the thermal broadening of the states at T temperature with ks the Boltzmann 
constant. 

Chen’s approach is based on the expansion of the tip wavefunction into spherical harmonic 
components around the tip apex position ro Q, 


Xv (r) = ^2 C uim ki («„r) Y lm (d, <p), (4) 

lm 

where r = |r — ro|, fa is the spherical modified Bessel function of the second kind, Yi m is the 
spherical harmonic function depending on the azimuthal ($) and polar (cp) angles and k v is 
the vacuum decay of the tip wavefunction. The expansion can also be performed according 
to the real space orbital characters (3 G {s,p y ,p z ,p x ,d xy ,d yz ,d 3z 2 _ r 2 ,d xz ,d x 2 _ y 2 } introducing 
the notation of Y v p{r, id, tp) = kp{^ l ,r)Yp{ , d,ip) with Yp real spherical harmonics as 


(r) = C »pY v p (r, p). 

0 


( 5 ) 


Using this expansion of the tip wavefunction in Eq. (j2J) leads to Chen’s derivative rule, and 
|M mv | 2 can be written as 


I M pv 


2 


47T 2 h 4 

K 2 m 2 


2 


0 ) 

P 


( 6 ) 


We introduce a new notation, M pv p = C v pd v p^) p { ro) that corresponds to the tunneling matrix 
element of a given orbital symmetry (/?). Here, the differential operator d u p acts on the sample 
wavefunction at the tip apex position tq. Note that d v p is dimensionless as it contains a factor 
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Table I. Differential operators d v p for given orbital symmetries (/?) according to Chen 


B- 


with l being the angular quantum number. The differential operators for the given orbital 
characters are summarized in Table [I] following Chen 7|. 

Rewriting Eq. (J6j) as 


= 


Ait 2 ?!* 
K 2 m 2 


0 0 ' 


Ait 2 ??* 

K 2 m 2 


Y I M-ni/0 1 2 + Y, 2 R e {M^ vj3 M^^} 

. 0 0 + 0 ' 


( 7 ) 


allows the investigation of different contributions to the tunneling current. The first term 
of the equation on the right hand side is the sum of the absolute value squares of M^ v p = 
Cu^dy^nir o), which is always positive, hence this term provides a positive contribution to the 
tunneling current. The second term is an interference term concerning tip orbitals, which is 
real and the sign of the individual (3^/3' components can be positive or negative, respectively 
contributing as constructive or destructive interference to the tunneling current. The analysis 
of the ratios and polarities of the listed components of \M^ U \ 2 gives the opportunity to obtain 
a deeper physical understanding of the electron tunneling process. Note that in Eq. ([6]) the 
sample wavefunction by can also be expanded into spherical harmonics similarly to Eq. (J5|). 
This way the interference of the sample orbitals and interference between sample and tip 
orbitals can be investigated as well. A similar decomposition of tunneling matrix elements 


has been used by Jurczyszyn and Stankiewicz 


20 


211 and Mingo et al. 4|. 


A. Calculation of spatial derivatives 


The spatial derivatives (see Tabled]) of the sample wavefunction for M^p = C^pd^p'ijj^r o) 
can be calculated straightforwardly when using a plane wave expansion of the wavefunctions 
There are many DFT codes, which use a plane wave basis set, e.g., VASP 


and Quantum-Espresso 


28j, ABINIT (291 


30] to name a few popular ones. Thus, the presented forms of the 


spatial derivatives can be potentially useful for future implementations of the revised Chen’s 
derivative rule. In the present work we use wavefunctions obtained from the VASP code. 
Let us assume that the single electron wavefunctions of the sample surface are given in the 
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vacuum at position vector r in a two-dimensional (2D) Fourier-grid as 


VVM =^nSkjf<rs( r ) = eXP [*( k f + G ||) r ||], 


( 8 ) 


where /! = (n^kj^tr 5 ) is the composite index for single electron states of the sample with 
kj| = (k^ky) the lateral component of the wave vector. The derivation with respect to z 
(the direction perpendicular to the sample surface) acts on the expansion coefficients only, 


d f d 

gZ^ n SkSas{ r ) = ( ^Z A nSk^(G||,^) ) exp [®(kjf + G||)r||], 


dz 


while the x- and ^-derivatives act on the phase factor, 
d 

;ts k ^s(r) = + G x )i4 nSk s <rS (G||,2)exp [i(kf H- G N )r N ], 


dx 

d 


tw( r ) = + G 'y)^k^(G||,^)exp [i(kjf + G N )r N ] 


dy^ nk \\ 

“ii 

The same procedure can be applied for higher order derivatives listed in Table Q] 


(9) 


( 10 ) 

( 11 ) 


B. Determination of weighting coefficients 


We report on three ways for the choice of the weighting coefficients C v p for M fW p = 

Cvpdvp^To). 

(i) The simplest choice is the assumption of an idealized tip with a given set of energy 
independent weighting factors {Cp}. Such examples can be found in the literature. In the 
study of Gross et al. a CO-functionalized tip was modeled as a combination of s and p orbitals 
and interference terms were neglected 14|. Siegert et al. employed the reduced density 


matrix formalism combined with Chen’s derivative rule with the inclusion of interference 

15|. Generally, Cp 


effects and they considered a similar combination of s and p orbitals 


can be a complex number. We restrict the choice of the set of { Cp } to fulfill the condition: 
Y2p \Cp\ 2 = 1. Moreover, in this idealized tip model case we choose the vacuum decay of the 
tip states k u =1 A -1 for all v. Examples of the effect of idealized tips as pure s and pure p z 
orbitals and a combination of (s + p z ) /\/2 on the STM image of N-doped graphene will be 
shown in section IIV Al We will also point out that the effect of interference is remarkable in 
this case causing a significant contrast change. 






(ii) Based on Eq. (J5]), C v p complex numbers can be obtained as 


C vP = 



{kp{Kur)Yp{d,(p) |xi/(r)) 


(12) 


with v = ( n 1 kjj 1 a T ) composite index for single electron states of the tip, where k| is the lateral 
component of the wave vector. We calculate these coefficients explicitly in the Wigner-Seitz 
sphere (W — S) of the tip apex atom with the VASP code. Since symmetry properties of 
the model tip geometry are taken into account in VASP, we obtain a reduced set of C u g 
corresponding to k| being in the irreducible part of the Brillouin zone. We can calculate 
how these coefficients change under 2D transformations (T) of the tip’s symmetry group in 
order to obtain C v p in the full 2D Brillouin zone. For this, the plane wave expansion of the 
tip wavefunction is needed, 


Xuir) = \v k [^(r) = ^B„T k T CT T(G||, 2 ;)exp [i(kjf + G||)r N ], (13) 

G u 

similarly to Eq. (J8]). Since the B expansion coefficients are invariant under the T transforma¬ 
tion, i.e., B n T k T a r = B T r T \ the transformation of the tip wavefunction comes from that 
ll n ' ( k n ) a 

of the phase factors. Using Eqs. (Eli and (TT3l) we obtain the following for the transformed 
coefficients, 

C nT T (k T )aTB = I kp(K v r)Yp(r) B V (G\\, z) exp [iT( kjf + G N )r N ] d 3 r 
V 111 M Jw-s n 

G ll 

= [ k0{K v r)Y0{Tr)xv{r)d 3 r. (14) 

Jw-s 

Note that T are represented by 2 x2 real matrices and the transformation of the coordinates is 
Tr = (7iiX+7i2?/, T 21 X+T 22 IIi z). Using the real spherical harmonics in Cartesian coordinates, 
we can calculate their transformations by substituting the transformed lateral coordinates 
into their normalized form. The results are shown in Table HH Thus, C v p is determined in 
the full 2D Brillouin zone, and we can directly apply them in the formula of the tunneling 
matrix elements in Eq. ()6l) . 

(iii) The third suggestion for C v p is based on the orbital-decomposed DOS projected 
to the tip apex atom, n TIP (E) = Yhp n ^ /P (-^) — Yhp ZL n up P fi(E ~ E v ) obtained from 
first principles calculation. Using the expansion of the tip wavefunction in Eq. (J5)) and the 
approximation of orthonormality for Y u 0 (r,d,(p) within the Wigner-Seitz {W — S) sphere of 


9 



Orbital 

Y(x,y,z ) 

Transformed orbital 

s 

i 

1\pK 

s 

Py 

1 flu 

2 \j 7T r 

TilPx + T22Py 

Pz 

1 j~3z 

2y 7r r 

Pz 

Px 

1 fix 

2 \j it r 

TllPx + Tl2Py 

dxy 
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1 /15 xz 

2 Y 7T 

'Tlldxz “I - Tl2dyz 

dx^—yZ 

1 fl5x 2 -y 2 

4 y 7r r 2 

(Ti\ — 7^) d x 2_ y 2 + 2TnT 12 4, 


Table II. Transformation of real spherical harmonics under 2D symmetry operations ( T ) of the tip. 


the tip apex atom, { Y u p(r, p>) Y u pi(r, $, ip) the following is obtained, 

/ w—s 


n 


TIP ( 


w = E E nT i F ^ E - B -) = E (*-i 

v p 1 / v p 

(15) 

Thus, we can approximate the complex C v p coefficients with real values, C v p ~ \J in ^p P • This 
way Eq. (JUJ) is recast to 

47T 2 h 4 ' 


\M^\ Z = 


n 2 m 2 




(16) 


Since the calculation of the orbital-decomposed atom-projected DOS is routinely available 
in DFT codes, the presented approximation applied to the tip apex atom gives a widely 
accessible choice for the weighting coefficients in the revised Chen’s derivative rule. In section 
IIV Al we will demonstrate that the STM images obtained by the C v p ~ \Jn PI ^ p approximation 
provide good agreement with those calculated using the proper complex C u p coefficients 
according to Eq. (Eli. 


C. Inclusion of arbitrary tip orientations 

Since the electronic structures of the sample surface and the tip are generally calculated 
independently to allow more flexibility with their geometries, arbitrary orientations of the tip 
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can be included into the revised Chen’s method. This can be done by redefining the spatial 
derivatives of the sample wavefunctions corresponding to the orbital characters in the rotated 
coordinate system of the tip with respect to the sample surface. This rotation is described 
by a coordinate transformation, which is represented by a 3 x 3 matrix 1Z with elements 7 Zj. 


We use the explicit form of 1Z as in Refs. 


26 


m, 


n = 


cos ipo cos 'bo — sin ipo sin bo cos i9 o cos <p o sin bo + sin <pg cos bo cos i9 o sin <p o sin i?o 

— sin (fo cos bo — cos ipo sin bo cos i?o — sin tpo sin bo + cos ipo cos bo cos $o cos tp o sin 
sin bo sin — cos bo sin i?q cos i9 o 


(17) 


with the Euler angles (i9o, bo, ipo)- Using the Einstein summation convention, the relationship 
between the two set of coordinates, the rotated tip coordinates G {x', y', z!} and the sample 
coordinates x l G {x,y,z}, is the following: 

c} T tj c) n r i 

x' j = —-X i = iz\x % \ x i = —x'l = (TZ- 1 ) . x' j . (18) 

dx l 1 dx'i K K ’ 


With the help of these, we can relate the derivatives of the sample wavefunction ip with respect 
to the rotated tip coordinates x /J to the derivatives with respect to the sample coordinates 
x l as 


dip dip dx l dip , iy 

dx ,J dx i dx'i dx i ' J 


(19) 


Similarly, the second derivatives are 

d 2 ip 

dx^dx'i 


( d 2 iP \ 
\dx l dx i ) 




( 20 ) 


Using Eqs. (1191) and (1201) the transformed & v ^ differential operators corresponding to the 
rotated tip coordinate system can be constructed and employed in Eq. (|6]) for the tunneling 
matrix elements. Since the transformation is linear, this, in turn, results in redefined C v p 
weighting coefficients in Eq. (J6|) for the d u p operators given in the coordinate system of the 
sample listed in Table 01 


III. COMPUTATIONAL DETAILS 

Using the revised Chen’s derivative rule implemented in the BSKAN code |4 3, STM 
imaging of two functionalized surfaces of current interest is investigated: N-doped graphene 
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and an antiferromagnetic Mn 2 H complex on the Ag(lll) surface in combination with several 
tip models. Geometrical relaxations and electronic structure calculations of the surface and 


28| employing the projector 


tip models were performed separately using the VASP code 
augmented wave (PAW) method 311. 

N-doped graphene is modeled as a free-standing single-layer graphene sheet in a 7 x 7 
surface unit cell following Ref. Q and 16 A-wide vacuum perpendicular to the surface to 
avoid unphysical interactions between neighboring slabs. One carbon atom is replaced by 
nitrogen in the given superccll. The generalized gradient approximation (GGA) and the 
exchange-correlation (XC) functional parametrized by Perdew and Wang (PW91) { 32 1 were 
used together with a plane wave basis set energy cut-off of 400 eV and an 11 x 11 x 1 
Monkhorst-Pack |33] k-point sampling of the Brillouin zone. We found a planar lattice 


structure after geometrical relaxation following N-doping, in agreement with Ref. 


191. 


For the revised Chen’s derivative rule, idealized model tips of pure s, pure p z and a 
combination of (s + p z )/V 2 orbitals are initially considered, see section III Bl (i) for details. 


Since N-doped graphene has been experimentally probed with tungsten tips 


19j, we consider 


three tungsten tip models with different sharpnesses and compositions: Wbi un t, W s w r and 
Wc-apex- The Wbiunt tip model is represented by an adatom adsorbed on the hollow site 
of the W(110) surface, the W s harp tip is modeled as a pyramid of three-atoms height on 
the W(110) surface, and the Wc-apex tip is a sharp tungsten tip with a carbon apex atom 
accounting for a likely carbon contamination from the sample. More details on the used tip 
geometries can be found in Ref. [10 |. 


Geometrical relaxations, search for the magnetic ground state and electronic structure 
calculations of an Mn monomer, Mn dimer and Mn 2 H on the Ag(lll) surface have been 


reported in Ref. 


22], and an antiferromagnetic ground state for the Mn 2 H/Ag(lll) system 


has been found. We use their electronic structure results in the present paper, for more 


details on the modeled geometries and DFT calculations please refer to Ref. 


22]. In their 


STM experiments silver tips have been used. Therefore, we consider blunt tips as an adatom 
adsorbed on the hollow site of the silver surface in two different orientations: Ag(001) and 
Ag(lll), in a 3 x 3 surface unit cell and at least 15 A-wide vacuum perpendicular to the 
surface to avoid unphysical interactions due to the slab geometry. GGA and XC functional 
parametrized by Perdew, Burke and Ernzerhof (PBE) 34} | are employed. Moreover, a plane 
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wave basis set energy cut-off of 250 eV and an 11 x 11 x 1 Monkhorst-Pack 


33] k-point grid 


centered on the T point are used. The convergence criterion for the forces acting on relaxed 
atoms (adatom and first full layer) is 0.01 eVA -1 . 


IV. RESULTS AND DISCUSSION 


We demonstrate the reliability of the revised Chen’s derivative rule for the mentioned 
N-doped graphene and antiferromagnetic Mn 2 H complex on the Ag(lll) surface, where 
quantum interference effects play an important role in the STM imaging process. This 
demonstration is done by qualitative and quantitative comparisons of simulated STM images 
with corresponding results obtained by Tersoff-Hamann and Bardeen tunneling methods. 
Quantitative comparison is facilitated by calculating Pearson product-moment correlation 
coefficients between the STM datasets |25l. |27|. Importantly, we find that the revised Chen’s 
model is 25 times faster than the Bardeen method concerning computational time taking the 
same tunneling channels, while maintaining good agreement. The effects of electronic struc¬ 
ture, orbital interference and spatial orientation of the tip on the STM images are highlighted. 
Since the detailed analysis of quantum interference effects and arbitrary tip orientations in 
STM junctions are presently highly demanding and enormously time-consuming using the 
Bardeen method, the implementation of the revised Chen’s derivative rule in the BSKAN 
code j^i, 31 is a very promising tool for more efficient STM simulations providing a deeper 
understanding of a wide variety of physical phenomena in STM junctions, e.g., quantum 
interference and tip geometry effects. 


A. N-doped graphene 


Experimental studies have shown that in the constant height STM images of N-doped 
graphene the tunneling current above the N atom is significantly lower than above the neigh¬ 
boring C atoms [19]. At first sight this seems to be in contradiction with the fact that 


the density of states of the N atom is larger than that of the neighboring C atoms close 
to the Fermi level. The current dip above the N atom has been explained by a destructive 
interference between the orbitals of the N and the nearest neighbor C atoms, a pure sample 


effect 


19]. Such a quantum interference effect is an ideal candidate to study with the revised 
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Chen’s method. 


We have calculated constant height STM images of the N-doped graphene surface using 
four different tunneling models: 3D-WKB, Tersoff-Hamann, revised Chen and Bardeen. The 
constant height STM simulations were performed at relatively small tip-sample distance (4 
A) at two selected bias voltages (±0.4 V) corresponding to the STM experiments by Telychko 


et al. 


19|. First, the 3D-WKB method |8j] has been used. This model takes into account 


the orbital characteristics and electronic structure of the sample and the tip as well, but 
uses the atom-projected density of states (amplitudes) instead of the explicit wavefunctions 
(amplitudes and phases), thus electron interference effects are not considered. Using the 
3D-WKB method the N atom always shows up as a protrusion in the STM image as it is 
expected from the relation of the density of states of the N and C atoms |l9|. 



Figure 1. Constant height STM images of N-doped graphene at 4 A tip-sample distance and ±0.4 
V bias. Comparison between Tersoff-Hamann model and revised Chen’s method with selected tip 
orbitals: s, p z , (s + p z )/V 2. The N defect is located in the middle of the images. 


Next, we focus on the comparison of the revised Chen’s method with two conventional 
STM simulation models: Tersoff-Hamann and Bardeen. Fig. |T] shows that the STM images 
obtained by the revised Chen’s method using a pure s tip quantitatively agree with those cal¬ 
culated by the Tersoff-Hamann model, i.e., a correlation value of 1 between the corresponding 
STM images is found. We obtain qualitatively similar STM images assuming an ideal tip of 
pure p z orbital, where the current dip above the N atom is slightly more pronounced than 
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with the s tip. For these tip models the energy independent weighting factors in section 
Kb] (i) were used. Furthermore, we point out the importance of quantum interference of 
tip orbitals in Fig. [TJ Therefore, we consider an ideal tip of a linear combination of s and 
p z orbitals of equal weights. As can be seen, the (s + p z )/V 2 tip shows a remarkable con¬ 
trast change compared to the STM images of pure s or pure p z orbitals, where the current 
dip above the N atom is even more pronounced and the bright triangle showing C atoms 
is larger and reversed. This effect is clearly related to the orbital interference of the tip’s s 
and p z states, and shows an additional effect to the destructive quantum interference arising 
from the sample’s C-N bond in the formation of the STM contrast of N-doped graphene. 
We stress again that the s — p z tip orbital interference results in a much more pronounced 
current dip above the N atom than the destructive quantum interference of the sample itself, 
the latter is imaged by the Tersoff-ffamann method. Interestingly, STM images obtained by 
the (s + p z )/V 2 tip resemble results calculated by a C(lll) tip model (see Ref. [lfij]) having 
these dominant orbitals in the electronic structure. Note that both types of STM contrast 

n 

of N-doped graphene in Fig. [l]have been experimentally observed in Ref. [19l |. 

Fig. [2] and Table UTT1 show qualitative and quantitative comparisons of the revised Chen’s 
method with Bardeen’s method using three different tungsten tip models, which have also 


been used in previous studies of STM imaging of HOPG 
different choices of the C v p weighting coefficients in Eq. 


10 


, |25]. Moreover, we compare two 
of the revised Chen’s method, 


see section KB] (ii) and (iii), and good qualitative agreement is obtained. 

We find that using these methods the current dip above the N atom is always present in 
the STM images in Fig. [2] The degree of agreement between Bardeen’s and revised Chen’s 
method, reported as correlation coefficients between corresponding STM images in Table HIT] 
depends on the actual tip geometry and electronic structure, hence the bias voltage. Let us 
recall that we expand the tip wavefunctions/density of states around the tip apex atom and 
calculate the C u p coefficients in the Wigner-Seitz sphere of the tip apex atom in the revised 
Chen’s method. The accuracy of such an expansion depends strongly on the neighboring sub¬ 
apex atoms’ electronic structure and on the tip-sample distance. For example, for sharp tips 


the con 
tips 


ribution of sub-apex atoms to the tunneling current is more important than for blunt 


25| . On the other hand, the larger the tip-sample distance the better the agreement 


of STM images between the two methods. The reason is that with increasing tip-sample 
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Figure 2. Constant height STM images of N-doped graphene at 4 A tip-sample distance and ±0.4 V 
bias. Comparison between Bardeen’s model and revised Chen’s method with two different choices of 
the C u p weighting coefficients in Eq. ()6]) (see section III Bl (ii) and (iii) for details) for three tungsten 
tip models: Wbiunt) W s } la rn and Wc-apex- Room temperature was assumed corresponding to the 


STM experiments in Ref. 


19y . The N defect is located in the middle of the images. 


separation the effect of the local tip geometry decreases. At larger tip-sample distances we 
find that the current dip above the N atom vanishes and a rounded triangular pattern is 
obtained leaving the three nearest neighbor C atoms visible, similarly to the Tersoff-Hamann 
results in Fig. [TJ 


Correlation coefficients 

Wbiunt tip 

Wsharp tip 

Wc-apex tip 

between STM images in Fig. [2] 

-0.4 V 

±0.4 V 

-0.4 V 

±0.4 V 

-0.4 V 

±0.4 V 

Bardeen-Revised Chen {C v p in Eq. (1121)1 

0.81 

0.82 

0.78 

0.89 

0.73 

0.91 

Bardeen-Revised Chen (C v p ~ \J n ^ P ) 

0.71 

0.81 

0.62 

0.79 

0.74 

0.92 


Table III. Quantitative comparison between Bardeen’s model and revised Chen’s method: calculated 
correlation coefficients between STM images in Fig. [2] 
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Overall, we find good agreement between Bardeen’s and revised Chen’s method in Fig. 
[2] and Table IIII1 Calculated correlation values are above 0.7 except for the W s harp tip at 
—0.4 V bias and C v p s= 


n 


TIP 
up • 


We find better correlation for C v p in Eq. (112D than for 
the C u q ~ approximation used in the revised Chen’s method, with the exception of 

the Wc-apex tip. Moreover, systematically better correlation between Bardeen’s and revised 
Chen’s method is found for +0.4 V than for —0.4 V bias voltage. Note that the large difference 
between STM images of the W s harp tip with different bias polarities can be explained by the 


asymmetric electronic structure of the tip apex around the Fermi level 
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Figure 3. Effect of strain on constant height STM images of N-doped graphene at 4 A tip-sample 
distance and ±0.4 V bias using the revised Chen’s method with three tungsten tip models: Wbiuntj 
W s h ar p and Wc-apex- The ground state N-doped graphene geometry obtained by DFT and two 
other structures with bond lengths varied by ±10% relative to the ground state are compared. 


It is interesting to investigate the effect of strain on the obtained STM contrast. Fig. [3] 
shows a comparison between STM images calculated for three different N-doped graphene 
geometries with varying bond lengths by ±10% relative to the ground state structure, which 
has been obtained by DFT calculation with C-N and C-C bond lengths of 1.42 A. Generally, 
we observe that the main features of the STM contrast do not change with the applied strain. 
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This is quantitatively confirmed by correlation coefficients being above 0.93 for each tip and 
bias combination calculated between images within each column of Fig. [3j We find a tendency 
of spatially extended brighter features in the STM images upon elongation of the bonds. 



Figure 4. Tip rotation effect on the constant height STM images of N-doped graphene calculated 
with W s harp tip at 4 A tip-sample distance and +0.4 V bias voltage. The rotation axis of the tip is 
perpendicular to the surface ($o = Ao = 0°). The orientations of the brightest features are indicated 
by white arrows in each STM image. 


In the following we investigate the effect of tip rotations on the STM images. Constant 
height STM images have been calculated above the N-doped graphene surface with the tung¬ 
sten tip models at 4 A tip-sample distance. We considered tip rotations around the axis 
perpendicular to the sample surface. This corresponds to $ 0 — 0°, and in this case rotations 
with respect to + 0 and i/jq are equivalent, thus, we fixed i/jq = 0° and rotated the tip by (p 0 in 
10° steps. Since a more asymmetric tip is expected to have a larger tip rotational effect, we 
present results obtained by the W s h a r P tip at +0.4 V bias voltage. Selected STM images are 
shown in Fig. [JJ 

Due to the symmetry of the tip, the same image is obtained for (p 0 = 180° as for 
+o = 0°. We find that the current dip above the N atom is always present independently of 
the degree of tip rotation by ip 0 , but the intensity of the current above the surrounding C 
atoms changes with the tip rotation. There are certain directions denoted by white arrows 
in Fig. U where the brightest features occur that correspond to the largest currents above or 
close to nearest neighbor C atoms. These indicated directions rotate twice faster than the tip 
rotation by + 0 itself. The finding that such kind of tip rotations, where the z axis of the tip 


is not tilted with respect to the z axis of the sample ($0 = 0°), affect the secondary 
of the STM image is in agreement with previous results using the 3D-WKB method 


ea 


25 


ures 


26|. 


In Fig. [5] we extracted line sections of the constant height STM images presented in Fig. 
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Figure 5. Current profiles along x and y directions of the N-doped graphene surface (see inset) as a 
function of tip rotation with respect to tpo- ($o = ipo = 0 °) 

[4j The line along the x direction contains the N atom and its nearest (Ci) and third nearest 
neighbor (C 3 ) carbon atoms. The line along the y direction contains the other two nearest 
neighbor carbon atoms (C) and C"), see the inset of Fig. [5^). The symmetries of the sample 
and the tip are reflected in these line sections as well. We find indeed that the current value 
above the N atom is insensitive to the tip rotation, and it is almost the smallest current value 
in the entire scan area. We can also see that the brightest features, i.e., the largest current 
values of the STM images are actually not located above the carbon atoms, but rather above 
the hollow positions of the honeycomb lattice. 


B. Mri2H on Ag(lll) surface 


Sachse et al. found that Mn 2 H on the Ag(lll) surface can produce STM images with 


single or double features depending on the magnetic coupling between Mn atoms 


22]. Double 


features have been obtained at positive bias employing the Tersoff-Hamann method for an 
antiferromagnetic Mn-Mn coupling, which corresponds to the energetically favored ground 
state. The calculated relaxed geometry of antiferromagnetic Mn 2 H on the Ag(lll) surface is 
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Figure 6. Calculated relaxed geometry of antiferromagnetic M 112 H on the Ag(lll) surface. Data is 


taken from Ref. 


22 ]. 


shown in Fig. [6] We consider this system and perform an investigation of its STM imaging 
depending on three employed tunneling models: Bardeen, revised Chen and Tersoff-Hamann. 
Using the decomposition of the tunneling current according to Eq. (J7J) in the revised Chen’s 
method, we are able to identify the physical origin of the observed dip above Mn 2 H. 


The calculated constant height STM images at small bias voltages (±0.1 V) using two 
silver blunt tip models [Ag(001) and Ag(lll)] are shown in Fig. 0 Correlation coefficients 
between STM images obtained by Bardeen’s and revised Chen’s method are reported in Table 
m First of all, we find excellent quantitative agreement between the STM images obtained 
by Bardeen’s and revised Chen’s method for the Ag(lll) tip and good agreement for the 
Ag(001) tip. Recalling that the revised Chen’s method is 25 times faster than Bardeen’s 
method in practical STM calculations, this clearly indicates that our proposed model is a 
very promising tool for STM simulations in the future. Moreover, the results in Fig. 0 show 
that the geometry and electronic structure of the tip have a considerable effect on the STM 
imaging of Mn 2 H/Ag(lll): Ag(001) tip provides single protrusion and Ag(lll) tip provides 
double features of the STM images at both bias polarities using both Bardeen’s and revised 
Chen’s methods. However, the Tersoff-Hamann model provides qualitative agreement with 
these at selected bias voltages only: At —0.1 V a single protrusion is obtained, while at 
±0.1 V a double feature is visible. The comparison of the Tersoff-Hamann results with those 
obtained by the revised Chen’s method suggests a contradiction with the general assumption 


of Ag tips being of s orbital character [22]. In order to understand the components of the 
tunneling current above the H atom, the decomposition according to Eq. (]7|) in the revised 
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Figure 7. Constant height STM images of antiferromagnetic M 112 H on the Ag(lll) surface simulated 
at 5 A Ag surface-tip distance and ±0.1 V bias with three different methods (Bardeen, revised 
Chen, Tersoff-Hamann) and two blunt tip models: Ag(001) and Ag(lll). A temperature of 7 K 


was assumed corresponding to the STM experiments in Ref. 


221. Note that results of the Tersoff- 


Hamann model are shown for comparison reasons only, and no explicit tip electronic structure is 
considered there. 


Chen’s method is employed. 

Fig. [8] shows the results of the current decomposition according to tip orbital characters. 
Interestingly, we find that the Ag(001) tip does not behave as an s-type tip at ±0.1 V bias 
[see Fig. [HR)]• The major contributions are from the p x , p y , d xz and d yz tip orbitals, and there 
are destructive interferences arising from p x — d xz and p y — d yz tip orbitals. This explains 
the qualitative disagreement between the Tersoff-Hamann and revised Chen’s results for the 
Ag(001) tip at ±0.1 V bias. On the other hand, using the Ag(lll) tip at —0.1 V bias, 
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Correlation coefficients 

Ag(001) tip 

Ag(lll) tip 

between STM images in Fig. [7] 

-0.1 V 

+0.1 V 

-0.1 V 

+0.1 V 

Bardeen-Revised Chen (C v p in Eg. (1121)1 

0.72 

0.88 

0.97 

0.95 


Table IV. Quantitative comparison between Bardeen’s model and revised Chen’s method: calculated 
correlation coefficients between STM images in Fig. [7] 


a) 


b) 



Figure 8. Decomposition of the tunneling current 1.83 A above the H atom in Mii 2 H/Ag(lll) 
(corresponding to Fig. 0 using Eq. ©• Diagonal: direct (positive) contributions, off-diagonal: 
interference (positive or negative) contributions to the current, a) Ag(001) tip at bias voltage 
V = +0.1 V, b) Ag(lll) tip at bias voltage V = —0.1 V. 

the major contribution is clearly from the tip’s s orbital [see Fig. [8b)]. Apart from that 
there is a strong s — p z destructive tip interference that is missing in the Tersoff-Hamann 
model, causing the observed qualitative difference in the STM images for the Ag(lll) tip 
at —0.1 V bias. Moreover, we find similar current decomposition characteristics for the 
Ag(001) tip at —0.1 V bias and for the Ag(lll) tip at +0.1 V bias as Fig. [8b) shows, with a 
dominating s orbital contribution from the tip. For these tip and bias voltage combinations a 
good qualitative agreement of the STM images between revised Chen’s and Tersoff-Hamann 
results is obtained. Our findings suggest that although the quality of the STM contrast 
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(single or double feature) is mainly determined by the electronic states of the sample surface 
that can be captured by employing the Tersoff-Hamann model, the tip electronic structure 
and in the present case an s — p z destructive tip interference can cause a contrast change. 



Figure 9. Effect of temperature on constant height STM images of antiferromagnetic M 112 H on the 
Ag(lll) surface simulated at 5 A Ag surface-tip distance and ±0.1 V bias using the revised Chen’s 
method with two blunt tip models: Ag(001) and Ag(lll). Temperatures of T = 7 K and 70 K are 
compared. 


It is important to highlight the effect of temperature on the obtained STM contrast. 
Temperature enters into the tunneling model exactly in the same fashion as in the well es¬ 
tablished Bardeen’s model, i.e., in two ways: (i) the thermal broadening of the electron states, 
see Eq.(j3J), and (ii) the energy window for calculating the tunneling channels, where nonzero 
temperature results in an extension of the energy window due to the Fermi distribution, see 
Eq.ffTlh Fig. [9] shows a comparison between STM images calculated at T = 7 K and 70 K, 


the former corresponding to the temperature used in the experiments of Ref. 


22). We find 


a diverse behavior of the STM contrast at the higher temperature depending on the tip and 
bias voltage. The contrast (single protrusion) is preserved for the Ag(001) tip at —0.1 V only. 
The other three images show different contrasts at the two temperatures. Upon increasing 
the temperature, for the Ag(001) tip at ±0.1 V the single protrusion contrast changes to 
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double features, while for the Ag(lll) tip at both bias voltages the double protrusion con¬ 
trast changes to an elongated single feature with the maximal current above the H atom of 
Mn 2 H. This diversity of simulated STM contrasts points to the importance of the correct 
choice of temperature in STM simulations if a meaningful explanation of given experimental 
STM data is desired. 


V. CONCLUSIONS 


We revised Chen’s derivative rule for electron tunneling for the purpose of computation¬ 
ally efficient STM simulations based on first principles electronic structure data. The revised 
Chen’s model includes the electronic structure and arbitrary spatial orientation of the tip by 
taking appropriate weighting coefficients of tunneling matrix elements of different tip orbital 
characters. Interference of tip orbitals in the STM junction is included in the model by 
construction. We demonstrated the reliability of the model by applying it to two function¬ 
alized surfaces of recent interest where quantum interference effects play an important role 
in the STM imaging process: N-doped graphene and an antiferromagnetic M 112 H complex 
on the Ag(lll) surface. We found that the revised Chen’s model is 25 times faster than the 
Bardeen method concerning computational time, while maintaining good agreement. Our 
results show that the electronic structure of the tip has a considerable effect on STM images, 
and the Tersoff-Hamann model does not always provide sufficient results in view of quantum 
interference effects. For both studied surfaces we highlighted the importance of interference 
between s and p z tip orbitals that can cause a significant contrast change in the STM images. 
Moreover, our findings show that stretched bonds have a minor effect on the main features of 
the STM contrast, and temperature is an important factor to be taken into account in STM 
simulations if aiming at accuracy in comparison with experiments. Our method implemented 
in the BSKAN code, thus, provides a fast and reliable tool for calculating STM images based 
on Chen’s derivative rule taking into account the electronic structure and local geometry of 
the tip apex. 
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